###############################################################################
#########################14. Robustness Droping each country###################
###############################################################################

#####Load data from Data_Generation_CPS#####

pols <- readRDS('pols_bjps')

#####Load plm package#####

library(plm)
library(pcse)
library(lmtest)
library(dplyr)
library(tidyr)
library(ggplot2)
library(reshape2)
library(rio)

#####Create Loop to capture T Stats of relevant variables#####

country <- as.character(unique(pols$country.name)) #List of countries

country <- country[-which(country == 'Cyprus' | country == 'Luxembourg' | country == 'Malta')] 

data <- matrix(nrow = length(country), ncol = 4) #matrix to hold output from loop

data[,1] <- country #fill in first colum of loop


for(i in 1:nrow(data)){
  m <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc1 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc1 + 
             scale_dispgini:scale_socioeconomicsal +
             scale_ggini_inc1: scale_socioeconomicsal +
             scale_ggini_inc1:scale_socioeconomicsal : scale_dispgini +
             scale_galtansdw +
             scale_enep +
             scale_unemp,
           data = pols[pols$country.name != data[i,1],], 
           na.action = na.omit, 
           model = "within")#Run model excluding country i
  
  #Save coefs adjusted for clustering#
  
  a<- as.data.frame(summary(m, vcov = function(x) vcovBK(x, cluster="group"))$coefficients)
  
  data[i, 2] <- a$`t-value`[1] #T-stat on gini main effect
  data[i, 3] <- a$`t-value`[7] #T-stat on gini*sorting
  data[i, 4] <- a$`t-value`[8] #T-stat on gini*salience
}

####Make plot of T-Statistics#####

plot_data <- data.frame(country = data[,1], 
                        disp_gini = as.numeric(data[,2]),
                        ginisort = as.numeric(data[,3]),
                        ginisal = as.numeric(data[,4])) #make data a data frame

melted_data <- melt(plot_data, id = "country")

melted_data$value <- as.numeric(melted_data$value)

#Plot data#

drop_test <- ggplot(melted_data, aes(x = country, y = value, shape = variable)) +
  geom_point() +
  theme_bw() +
  scale_y_continuous(limit = c(0,7)) +
  labs(x = "Country Dropped from Model", 
       y = "T - Statistics") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom") +
  scale_shape_discrete(name = "Variable:", labels = c("Disposable Income Inequality", "Interaction of Inequality and Partisan Sorting", "Interaction of Inequality and Economic Salience")) +
  geom_hline(yintercept = 2, linetype = 'dashed')

#ggsave("figure_A2.pdf", plot = drop_test, width = 11, height = 9, units = "in" )
